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ABSTRACT 



Context. With the current Herschel and Planck satellite missions, there is strong interest in the interpretation of the details of the 
sub-millimetre dust emission spectra from interstellar clouds. A lot of work has been done to understand the negative correlation 
observed between the spectral index /?obs and the colour temperature Tq that in the;^^ fits is partly caused by the observational noise. 
Aims. In the (Tq, jSqu) plane, the confidence regions are elongated, banana- shaped structures. Previous studies have indicated that 
the errors may exhibit strongly asymmetric features that have important consequences for the investigation of individual objects and 
the interpretation of the relation between the Tq and/?obs parameters. We study under which conditions the confidence regions exhibit 
such anomalous, strongly non-Gaussian behaviour that could aff'ect the interpretation of the observed (rc,/?obs) relations. 
Methods. We examined a set of modified black body spectra and spectra calculated from radiative transfer models of filamentary 
interstellar clouds. We analysed simulated observations at discrete wavelengths between lOOyt/m and 850yL/m. We performed modified 
black body fits and examined the structure of the ;t'^(rc,/?obs) function of the fits. 

Results. We demonstrate cases where, when the signal-to-noise ratio is low, the;^^ has multiple local minima in the (rc,/?obs) plane. 
A small change in the weighting of the data points can cause the solution to jump to completely diff'erent values. In particular, if there 
is noise, the analysis of spectra with T > lOK and/?obs ~ 2 can lead to a separate population of solutions with much lower colour 
temperature and higher spectral indices. The anomalies are caused by the noise. However, the tendency to show multiple ;^^ minima 
depends on the model (in part via the influence on the signal-to-noise ratios) and on the set of wavelengths included in the analysis. 
Deviations from the underlying assumption of a single modified black body spectrum are not significant. 

Conclusions. The presence of several local minima implies that the results obtained from the;^^ minimisation depend on the starting 
point of the optimisation and may correspond to non-global minima. Because of the strongly non-Gaussian nature of the errors, 
the obtained (Tc,j3ou) distribution may be contaminated by a few solutions with unrealistically low colour temperatures and high 
spectral indices. Proper weighting must be applied to avoid the determination of thcy^obsCT'c) relation to be unduly aff'ected by these 
measurements. 

Key words. ISM: clouds - Infrared: ISM - Radiative transfer - Submillimeter: ISM 



1. Introduction 

Observations of the thermal dust emission is one of the main 
methods used to map dense interstellar clouds. Sub-millimetre 
and millimetre emission is considered one of the most reliable 
ways of obtaining information about cloud cores in the dif- 
ferent stages of the star form ation process, from starless cores 
to protostellar sy stems (.Motte et al.l 119981; lAndre et al] I2QQQI: 
lEnoch et aP l2007h . To a large extent, this conclusion is based 
on the problems with the other tracers, i.e., the difficulty of in- 
terpreting molecular line data of the very dense clouds and the 
difficulty of assembling hi gh-resoliition maps using dust extinc- 
tion or scattering (.Lornbardi et al.ll2QQ6[ iGoodman et al.l [20091: 
lJuvelaetal.ll2QQ8[l20Q9L 

The main difficulties in the interpretation of dust emission 
are well known. Because of the line-of- sight temperature vari- 
ations, the peak of the emission spectrum becomes wider and 
this results in a decrease of the observed spectral index /5obs 
(IShettv et al . 2009b; Malinen et al. 201 1; Juvela & Ysard 201 1). 
At the same time, because the observed emission is dominated 
by the warmer dust components within the beam, the colour tem- 



perature Tq overestimates the mass averaged physical dust tem- 
perature. T his can lead to a s i gnificant underestimation of the 
dust mass (lEyans et al.ll200ll:lstamatellos & Whitworthll20Q3l: 



iMalinen et al.1 1201 ll: lYsard et al.1 1201 ll) . Furthermore, it is dif- 
ficult to estimate the intrinsic dust grain properties, such as the 
spectral index, on the basis of the observed radiation. 

When far-infrared and sub-millimetre observations are fitted 
with modified black body spectra, the dust colour temperature, 
Tq, and the dust spectral index, y^obs, are partially degenerate. A 
small increase in Tq can be compensated by a small decrease of 
y^obs- When the observations contain noise, the (rcy^obs) values 
will scatter over an elongated ellipse with a negative correlation 
between the two parameters. If the noise is large enough, the 
points will fall on a banana- shaped region where the Poh^jTc) 
relation becomes steeper at lower values of Tc ( Shettv et al. 
| 2009bllal: IVenezianietan 120101: iParadis et al.1 120101: 1 Juvela et al ' 



zOTT). The common assumption is that apart from the curvature 
of the error banana, the scatter is symmetric with respect to the 
true values of Tq and y^obs so that the expectation values do not 
show significant bias. 
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TheySobs(^c) relation induced by the noise is usually steeper 
than the average ySobsCT'c) relation derived from observations 
of a large numb e r of individual objects (P upae et al. 2003; 
iDesert et alJbOOSl: iPlanck Collaboration et al.l l201 1). If the av- 
erage of the observed (Tq, y^obs) points remains on the in- 
trinsic I3(T) relation and if one observes a wide range of ob- 
jects with different true temperatures, the effect of the bias 
can be controlled. However, there are some indications that 
under certain conditions, the error distribution is not symmet- 
ric and could behave in an even more non-Gaussia n fash- 
ion (iPlanck Collaboration et"aDl201 ll: lYsard et al.ll201 ih . Under 
some conditions, thQ/3ohs(Tc) can extend to low colour temper- 
atures and very high values of the spectral index. This could be 
important for the interpretation of the observed ySobs(^c) rela- 
tions and, more generally, any attempt to measure cloud masses 
and temperatures using noisy continuum data. 

In this paper we investigate this question by analysing noisy 
modified black body spectra, mixes of these spectra with dif- 
ferent colour temperature, and spectra obtained from radiative 
transfer modelling of optically thick clouds. The content of the 
paper is the following. In Sect. |2] we describe the methods used 
to produce the spectra and to derive the colour temperature and 
the spectral index estimates. The main results are presented in 
Sect. O We start by looking at spectra calculated for optically 
thick filaments and by identifying the cases where the x^ has 
more than one local minimum. In the Sect. l3.2l we return to mod- 
ified black body spectra in an effort to identify the basic require- 
ments for these anomalies. In Sect. |4] we discuss our results and 
the final conclusions are presented in Sect. |5] 



2. Methods 

We describe below the procedures used to calculate synthetic 
spectra and to derive the Tq andySobs values through ;^^ minimi- 
sation. 



2.1. Spectra from radiative transfer models 

The first set of model spectra was obtained from the radia- 
tive transfer models discussed in Ysard et al. (zOlT). The mod- 
els consist of optically thick filaments that are represented by 
long cylinders with radial density d istributions following the 
Tlummer-like' profiles ( Nutter et aDl20Q8l lArzoumanian et al.l 
l2Qllh . which are flat at the centre of the filaments and decrease 
as R~^ in the outer regions (R is the radius of the clouds). The 
central extinction is varied in the range Ay =1-20"^ and the 
analysed sub-millimetre emission remains optically thin. The 
clouds are externall y heated by the standard radiation field, ISRF 
(iMathis et al.ll 19831) . Because cloud cores are often embedded in 
large molecular cloud complexes, this radiation field can also 
be attenuated corresponding to an external layer of dust with 
Ay ^ = 1-5"^. We used dust properties representative of the dust in 
diffuse high Galactic latitude medium (DHGL) as defined in the 
DustEM dust models (Compiegne et al. 2011). The dust model 
consists of three dust populations: interstellar PAHs, amorphous 
carbons, and amorphous silicates. In the model clouds, the dust 
temperatures range from over 20 K for the amorphous carbon 
at the cloud surface to less than 10 K for the silicate grains at 
the centre of the filament, also depending on the model optical 
depth. Because of the wide range of temperatures and the dif- 
ferent intrinsic emissivity spectral indices JS of the dust popula- 



tionfl the spectra cannot be fitted precisely with a single mod- 
ified black body. For details of the calculations see lYsard et aP 
( 2011,) . 



2.2. Modified black body spectra 

We additionally examined spectra that are based on modified 
black bodies to which observational noise is added. The differ- 
ent cases are (1) a single modified black body with a fixed value 
of /3, (2) the sum of two modified black bodies with the same j3 
but different temperatures, and (3) a modified black body with 
different yS values below and above the wavelength of 500 yum. In 
the first scenario we are testing if the noise alone can produce 
a tail of solutions with high y^obs and low Tq values. With the 
other modifications, we tested if the probability of extreme val- 
ues is enhanced when the original spectrum cannot be described 
exactly as a single modified black body. 



2.3. Analysis of the simulated spectra 

In accordance with recent observational studies, we used mea- 
surements at a few far-infrared and sub-millimetre wave- 
lengths. To simulate the Planck studies, we used the wave- 
lengths of 350, 550, and 850 yum complemented with the 
lOOyum p oint that would be available fr om the IRAS sur- 
vey (e.g., iPlanck Collaboration et al.l |201 1]) . As a default, we 
added to the spectra observational noise that is 0.06, 0.12, 0.12, 
and 0.08MJysr-^ at the wavelengths of 100, 350, 550, and 
850 yum, respectively. For Planck the uncertainties of the abso- 
lute surface brightness measurements are significantly smaller 
but these numbers are more realistic if the flux determina- 
tion includes the separation of a background component (see 
iPlanck Collaboration et al.ll201 1\) . To simulate Herschel obser- 
vations, we used the wavelengths of 100, 160, 25 0, 350, and 
500 yum (e.g.. iParadis et aLlbOlCt iJuvela et al.ll201 ih with noise 
levels of 8.1, 3.7, 1.2, 0.85, and 0.35MJy sr'^ per beam. In the 
analysis the data are convolved to the resolution of the 500 yum 
observations and this results in final uncertainties of 1.62, 1.18, 
0.60, 0.60, and 0.35 MJy sr"\ respectively. After the beam con- 
volution, the absolute signal is lower for the simulated Planck 
data (a resolution of ~ 5') compared to the simulated Herschel 
data (a resolution of ~ 36'' at the wavelength of 500 yum). This 
reduces the difference in the signal-to-noise ratios (S/N) of the 
two data sets (see Fig.[T]). 

A weighted least- squares fit of a single modified black body 
(X^ minimisation) was used to estimate the colour temperature 
Tq and the observed spectral index /^obs- The wavelength ranges 
fitted are 100-8 5 Oyum for the combined data of Planck and IRAS 
and 1 00-500 yum for the simulated Herschel data. The weighting 
was done with the actual noise in each channel although the ef- 
fect of changing the weight of the lOOyum measurement was also 
examined. When there are temperature variations, the derived Tq 
andy^obs values differ from the mass weighted average dust tem- 
perature and from the actual spectral index of the dust grains. 
However, in this paper we concentrate on the effects of noise. 
Therefore we concentrated on the question of how the Tq and 
y^obs values differ from the values that would be obtained with a 
similar analysis of the noiseless spectra. 



^ The intrinsic spectral index of the amorphous carbon is 1.55, while 
it is equal to 2.11 for amorphous silicates. 
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Fig. 1. Signal-to-noise ratios in the simulated observations. The 
solid lines show the average and the dashed lines the minimum 
and the maximum S/N ratio as the function of wavelength. 




Fig. 2. Distribution of (T^c, y^obs) values for the modified black 
body fits of the spectra calculated for the filament models. The 
colour scale shows the density of points per Arc=0.2K and 
Ay^obs^O-l- The white circles (with black borders) indicate the 
values obtained for the examined models in the absence of noise. 



3. Results 

3.1. Spectra from radiative transfer models 

3.1 .1 . Simulated Planck and IRAS observations 

We start w i th a s tudy of the spectra that were calculated in 
lYsard et aP (12011") for models of cloud filaments at a distance 
of J=100pc. Figure [2] shows the (Tq, jSohs) values derived from 
simulated observations of 100, 350, 550, and 850yum with the 
default noise (see Sect 12.11) . In the figure we include models 
with Ay =10-20"^ and with the ISRF attenuated by an external 
dust layers with Ay^=4-5^. For these models the locus of the 
correct Tq andy^obs values (i.e., the parameters estimated in the 
absence of noise) is between the Tq values of 12.0 K and 12.8 K 
and the /3 values 1.02 and 1.11. Because of the noise the esti- 
mates scatter along a narrow banana-shaped region. Most points 
cluster around the median value of 12.47 K and yS= 1.07. The fig- 
ure contains 10000 points of which 217 are below Tc=^ K. 

The histograms in Fig. [3] show that the parameter distribu- 
tions are not symmetric and this is not caused merely by the 
curvature of the confidence region. The colour temperature dis- 
tribution has a tail towards low values and a secondary peak is 
visible around 7-8 K. The spectral index distribution is corre- 
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Fig. 3. Marginal distributions (i.e. the probabilities integrated 
over the j3 and Tq axes, respectively) of the points in Fig.lH The 
upper frames show the probability density distributions of the Tq 
and j3 parameters and the lower frames show a zoomed version 
of the upper frames. The histograms have been normalised to 
represent probability distributions (area normalised to one). The 
vertical dashed lines indicate the median of the distribution and 
the points where the tails of the distribution (one-sided) contain 
5% or 1% of the data. 



spondingly skewed in the other direction with a tail extending to 
high values ofjS. Altogether 4.6% of the points have Tq < 10 K. 
These represent a strongly non-Gaussian part of the error dis- 
tribution that, if not properly accounted for, would negate any 
attempts to determine what the realySobs(^c) dependence would 
be in the absence of noise. 

For filament models of lower Ay and thus of lower signal-to- 
noise ratio the high /^obs solutions become more common. The 
Ay = 20"^ models are responsible for -17% of all the points 
below Tc = 8 K, while the contribution of the clouds with a cen- 
tral Ay of 10"^ is already over 26%. More interestingly, the long 
tail to low colour temperatures is almost exclusively produced 
by the models where the external radiation field is attenuated by 
a dust layer of Ay ^ = 5"^. The models with Ay ^ = 4^ still show 
an asymmetry of the Tq distribution but their contribution to the 
tail below Tq = 10.0 K is only a couple of percent. Figure [2] 
also shows the (T^c, jSots) values that would have been obtained 
from the various models if there were no observational noise. 
The Ay ^ = 4"^ and Ay ^ = 5"^ models form separate groups, the 
latter being colder, as measured by the noiseless colour temper- 
ature, by AT^c ~1 K. The higher Ay ^ values reduce the physical 
dust temperature especially at the cloud surface. The eff'ect is felt 
most strongly at lOOyum where the signal-to-noise ratio drops by 
-40% (from ~6 down to ~3.5, almost irrespective of the central 
Ay of the model cloud). 

Figure |4] shows one example of a Ay = 5"^ cloud 
where the surface brightness values are 0.088, 1.62, 1.27, and 
0.35MJysr-i at the wavelengths 100, 350, 550, and 850yum, 
respectively. A fit to the original noiseless data gives values 
Tq =12.6 K and ySobs= 1.14. With this particular noise realisa- 
tion the;^^ minimum has moved to rc=4.78K andy^obs =4.46. 
These values are suspicious because the colour temperature is 
well below the minimum dust temperature of the model cloud, 
which is always higher than 8 K. The spectral index y^obs is also 
higher than the dust intrinsic j3 (1.55 and 2. 1 1 for amorphous car- 
bons and silicates, respectively, see Fig. 8 in lYsard et al.ll201 ih 
although it would be expected to be lower because of the line-of- 
sight temperature variations. In the model clouds the actual dust 
temperature does not decrease below ~8 K, while the spectral in- 
dices of the dust grains only are 2.11 for astronomi cal silicates 
and 1.55 for amorphous carbon grains (see Fig. 8 in lYsard et al.l 
■2011.) . 
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Fig. 4. Case where noise has resulted in a high /3 value. In the 
middle frame, the black curve is the spectrum obtained from the 
cloud modelling. The red symbols are the measurements that in- 
clude noise and the red curve is the fit to these points. The upper- 
most frame shows the deviations A/y (the diff'erence between the 
simulated observation and the fit) in units of the assumed uncer- 
tainty criy. The bottom frame shows the;^^ values as a function 
of Tc and/^obs- The colour plot and the white contours show the 
x'^ values for the fit to the noisy data. The contours are drawn 
at 1.02, 1.05, 1.1, 1.5, 2.0, 4.0, and 8.0 times the minimum ;^2 
value. In yellow are shown the corresponding contours for the fit 
to the noiseless data. The two circles denote the locations of the 
X^ minima with (red circle) and without (green circle) the noise. 



In this case the lOOyum intensity was not much more than 
l-cr detection. However, the lOOyum value happens to be almost 
identical to the correct noiseless value. On the other hand, the 
550yL/m observation is almost 2.5cr above the correct value. This 
is the main reason for the very low temperature. The contour plot 
of the ;^^ values (lower frame in Fig.|4l) shows that the confidence 
region is very elongated. The 'correct' solution (i.e., the one for 
the noiseless spectrum) resides in the x^ valley but with a x'^ 
value that is four times the;^^ value of the best fit. 

If the observations do not perfectly fit a single modified black 
body, the best fit will depend on the relative weight given to 
the individual measurements. Figure |5] shows that interesting 
changes take place when the assumed uncertainty of the lOOyum 
point is decreased by a factor of -0.57. The measured values 
(including the noise) are not changed and only the weight of the 
lOOyum point is increased in the fit. In the present example the 
lOOyum value happens to be almost correct. Therefore, the er- 
ror estimates could be decreased much more than by a factor of 
0.57 without this particular noise realisation becoming improba- 



ble. In Fig. [5] the lOOyum noise is scaled by 0.575 in the left hand 
frames and by 0.570 in the right hand frames. As the weight of 
the lOOyum point is modified, the;^^ minimum jumps instanta- 
neously from the high JS solution to a solution near the value of 
the original noiseless spectrum. At this point the signal-to-noise 
ratio of the 350 yum data point is ~18. The;^^ values are almost 
identical, the change being consistent with the eff'ect of the very 
small decrease of the lOOyum uncertainty. The parameter space 
was sampled with intervals Ar=0.1 K and AyS=0.02. 

The example reveals some important facts. Firstly, the so- 
lution based on x^ minimum can change significantly and in a 
non-continuous fashion when the measured intensities or the as- 
sumed uncertainties are slightly perturbed. Secondly, the pres- 
ence of multiple local x^ minima implies that the solution ob- 
tained by non-linear optimisation will depend on the initial val- 
ues of the optimisation. The obtained result may correspond to a 
local instead of a global minimum. Thirdly, as already indicated 
in Fig. O the error distributions are skewed and possibly even 
bimodal. This has implications not only for the uncertainties of 
the;^^ approach but more generally for the statistical models em- 
ployed in other parameter estimation methods. 

In the following we call 'normal' the solution that is close 
to the values obtained in the absence of noise and 'anomalous' 
the solutions that exhibit a significantly lower value of Tq and a 
higher value of y^obs- Another example is shown in Appendix lAl 
where the 550 yum point is more than 2cr above the correct (noise- 
less) value. In the 217 cases of the Ay > 10"^ model spectra with 
colour temperature below 8 K (-2% of all spectra), the common 
feature is that the 550 yum point is high relative to the neighbour- 
ing wavelengths and especially relative to the 850 yum point. This 
is demonstrated in Fig. [6] which shows the diff'erences between 
the observed and the true intensities when the estimated Tq was 
below 8K. With a low weight of the lOOyum measurement, a 
solution with a very high value of y^obs becomes possible. The 
anomalous solutions are seen more frequently but not exclu- 
sively when the lOOyum intensity is underestimated. One must 
also note that the previous colour temperature estimates (e.g.. 
Fig. [21) were based on x^ optimisation where the initial values, 
r =15 K andyS=1.5, favoured the high- temperature solution. 

It is time-consuming to estimate the x'^ values for the mod- 
ified black body fits with a dense grid over the whole (Tq, JS) 
plane. Instead, to identify the cases with two x^ minima, we 
started non-linear optimisation (the Powell method) at two loca- 
tions, (rc,yS)=(7.0K, 4.0) and (13.0 K, 0.9). If two minima exist, 
the optimisation is likely (although not guaranteed) to converge 
to diff'erent values. With these two initial values the optimisa- 
tion may still converge to the same local minimum, missing the 
second one. On the other hand, if there is only one very shal- 
low minimum, the optimisations may produce diff'erent results 
for numerical reasons. 

The presence of a single x^ minimum does not tell us 
whether it corresponds to the normal or the anomalous solu- 
tion. To see whether one is close to a situation where two x^ 
minima appear, we also scaled the lOOyum error estimates by 
a factor ^loo that was changed from 0.5 to 1.5 in steps of 0.1. 
This way one can perturb the problem and obtain an indication 
whether the result of the fit is well-defined or not. As seen in 
Fig. \5\ the minimum can shift very rapidly between the normal 
and the anomalous solution. The second local x^ minimum ex- 
ists only for few a ^loo values close to that point. Our emphasis 
is on the non-Gaussian nature of the uncertainties (as produced 
by the multiple minima). Note that for example Fig. [2] was ob- 
tained using only the original weighting of the data , ^loo =1.0. 
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Fig. 5. As Fig.|4]but assuming in the fit a 100 yum uncertainty that is 0.575 times (frames on the left) or 0.570 times (frames on the 
right) the original value. The intensity values are the same as before. The x^ plane exhibits two distinct local minima and a small 
change in the weight of the lOOyum measurement moves the solution from a high ft solution to a low j3 solution. 



We examined the above set of 10000 spectra from cylindrical 
cloud models with Av= 10-20"^ and A^f of 4-5"^. This revealed 
-680 cases where diff'erent initial parameter estimates lead to 
diff'erent optimised values with AT > 0.2 or A/5 > 0.1. This 
happened preferentially when the 100/im uncertainties were in- 
creased but could take place for any usually small range of ^loo 
values. In these models the ratio of the original intensity and 
the added noise was at lOOyum higher than 3.2, with an aver- 
age value of 4.7. Of the three Planck channels the 500 yum band 
had the lowest ratio with a minimum of 15.1 and a mean value 
of 18.7. The actual signal-to-noise ratios can be lower when the 
observed intensity is below the expectation value of the intensity. 

The double x^ minima are caused by the noise but their fre- 
quency of appearance depends on the model, probably mainly 
through the associated changes in the signal-to-noise ratios of 
the observations. For the same model clouds with Ay =10-20"^ 
but with the external radiation field attenuated by Ay ^ < 4"^ 
rather than 4-5"^, the double minima become less frequent by a 
factor often. With Ay ^ < 2^ the solutions were unique, probably 
because of the higher surface brightness values of those models. 



3.1.2. Simulated Herschel observations 

The simulated Herschel observations consisted of the wave- 
lengths of 100, 160, 250, 350, and 500 urn. Figure |7] s hows the 
(Tc,fi) values for the models from lYsard et al.l (1201 ih with the 
cloud central Ay =10-20"^ and an external Ay ^=4-5. T his is only 
a subset of all models discussed in lYsard et al.l (l201ll) . The orig- 
inal weighting of the lOOyum data has not been altered (i.e.. 
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A [ftm) 

Fig. 6. Errors in the observed intensities in the cases where the 
simulated Planck and IRAS lOOyum observations result in colour 
temperatures Tq <8 K. The errors are given as the diff'erence be- 
tween the observed intensity and the noiseless spectrum, mea- 
sured in units of the uncertainty assumed in the modified black 
body fits. 



^100 =1.0). The figure includes 10000 points out of which only 
83 are below Tc =9 K. 

The scatter of the temperature and spectral index values is 
more symmetric than for the Planck -I- IRAS data set. This is con- 
firmed in Fig. [8] which shows the marginal distributions where 
the Tq data have only a hint of a tail towards higher tempera- 
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Fig. 7. Distribution of (Tq^JS) for simulated Herschel obser- 
vations. The colour scale gives the density of the points per 
Arc=0.2K and A/5obs=0-l The white circles (with black bor- 
ders) indicate the values that would be obtained for the included 
models if there were no noise. 
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Fig. 8. Probability distributions of Tq and JS for the simulated 
Herschel observations of Fig. [7] (for further details, see the cap- 
tion of Fig. [3]). 



tures (skewness 0.63), i.e., opposite to the behaviour in Fig.O 
All spectra were again fitted using an optimisation with diff'erent 
initial values of Tq andyS. Out of the 10000 spectra with Av=10- 
20"^ and Ay H-5"^, two;^^ minima were inferred only in a single 
case. This even though the lOOyum S/N ratios was below one 
and the 160yum ratios only a few, the mean value being 5.0. For 
the longer wavelengths, the S/N ratios were ~20 or above. The 
X^ plane was examined for some of the cases with the lowest 
colour temperature to confirm the presence only of a single local 
minimum. The other two samples, the case with Ay =10-20^ and 
^ext ^ 4m ^^^ ^^Q ^^gg ^j^l^ Av=l^ and Ay ^ < 1"^, together con- 
tain only one additional case where two local x'^ minima were 
detected. 

On the basis of the previous tests the appearance of anoma- 
lous solutions depends mainly on the noise. The comparison of 
the Planck and Herschel cases shows that the set of wavelengths 
included in the analysis is equally important. The uncertainties 
of the simulated Herschel observations are all higher in absolute 
terms but the relative uncertainty between the short and the long 
wavelengths is similar to the Planck case. The combination of 
five wavelengths appears to be resistant against extreme errors 
that could be caused, for example, by a 2-cr or 3-cr error in a 
single channel (cf. Fig|4]). The Tq andy^obs distributions are wide 
because of the noise but there still are practically no values of 
Tc < 8 K. 

There is still the possibility that the results could depend on 
how much the underlying spectrum deviates from a single mod- 
ified black body. To further examine this question, we next ex- 
amined a series of models based on modified black bodies. 
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Fig. 9. Frequency of recognised;^^ double minima for simulated 
IRAS and Planck observations. Each frame corresponds to one 
combination of the temperature and the spectral index that were 
inputs of the simulation. The horizontal axis gives the 350yum 
intensity (before noise is added) and the vertical axis is the fac- 
tor ^100 (smaller numbers correspond to a larger weight of the 
lOOyum point). The colour scale gives the frequency of the dou- 
ble x^ minima as a percentage of all noise realisations with the 
corresponding values of S (350 mum) and ^100 • 



3.2. Modified biack body spectra 

We examined a series modified black body spectra By(T) x v^ to- 
gether with noise to look for similar anomalous cases as shown 
in Fig. [5] Several combinations of the temperature and spectral 
index were investigated. We started with the lOOyum band and 
the three Planck bands 350-850yum. When the modified black 
body spectra are scaled to have a 350yum signal of 2.0 MJy sr"\ 
the signal-to-noise ratio is comparable to that of Fig. [5] However, 
to examine the eff'ect of diff'erent S/N ratios, the spectra were 
scaled by a factor that was varied from 0.6 to 6.0 in five logarith- 
mic steps. To investigate the sensitivity to the relative weighting 
of the frequency points in the fit, we included the factor ^100 
which was varied from 0.5 to 1.5. As before, the factor ^100 did 
not aff'ect the noise, only the relative weighting of the data points 
in the fit. For each case (i.e., a combination of temperature, spec- 
tral index, the intensity scaling, and the value of ^100), we ran 
5000 noise realisations of the spectrum. Each spectrum was fit- 
ted using the two diff'erent initial values of the optimisation to 
recognise the cases with separate x^ minima. 

Double minima are detected in -10% of the cases but the fre- 
quency drops close to zero by the time the S/N ratio is increased 
by a factor of three (Fig. O. There is some dependence on the 
model in question but this may also be caused by the diff'erences 
in the S/N ratios of the other bands. With T = 8.0K andyS =1.5, 
the double minima are significantly more rare when the lOOyum 
data point is given a larger weight (^100 <1.0). With T = 12.0 K 
andyS =2.5, the opposite is true. 

Figure [TOl shows the corresponding results for the Herschel 
data that were scaled to have the same S/N ratios at the 350 yum 
wavelength as in the previous Planck examples (the observa- 
tional uncertainties are the same as in Sect. 13.11 ). Distinct x^ 
minima are observed only with Tq = 8.0 K andyS =1.5 and 
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Fig. 10. Frequency of recognised x'^ double minima for sim- 
ulated Herschel observations. Each frame corresponds to one 
modified black body with the given temperature and the spectral 
index. The horizontal axis gives the 350 yum intensity (the S/N 
ratios are the same as in Fig© and the vertical axis is ^loo, the 
relative weighting of the lOOyum point. The colour scale gives 
the frequency of the double ;^^ minima as a percentage. 
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Fig. 12. Frequency of the recognised double x^ minima for the 
sum of two modified black bodies with temperatures 8 and 10 K 
or 10 and 14 K. The column density ratio of the lower temper- 
ature and the higher temperature component is 1:2 or 2:1, as 
indicated in the frames. The data correspond to simulated obser- 
vations at lOOyum and three Planck wavelengths. 
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Fig. 11. Frequency of the recognised ;^^ double minima for sim- 
ulated Herschel observations when the signal-to-noise ratio has 
been decreased by a factor of seven compared to Fig. [TOl 



even in that case the number is below 1%. When the signal-to- 
noise ratio is increased by a factor of three, no double minima 
are detected (probability ~ 10"^ or less). 

If the 350yum signal of the simulated Herschel observations 
is lowered to 2.0MJysr"^ the S/N ratios decrease by a factor 
of seven. As shown in Fig.[TTl in this case the number of double 
minima increases to -10%, a level similar to the previous Planck 
example. 

The figures in Appendix |B] show the bias and scatter of the 
Tc and/^obs values corresponding to the cases in Figs.|9]-[TT] 



We next examined spectra that are the sums of two modified 
black bodies that have the same spectral index but diff'erent tem- 
peratures, /y = (By(Ti) -\- By(T2)) X v^ . Tfic deviations from a 
single modified black body could make the fit more susceptible 
to ambiguity. 

The results for the Planck-hlRAS case are presented in 
Fig. [121 They look quite similar to the single black body cases 
shown in Fig. [9] where, of course, also the temperatures are 
somewhat diff'erent. It seems that the deviations from a single 
modified black body shape do not have a significant eff'ect on the 
presence of multiple ;^^ minima. The conclusion is confirmed by 
the Herschel simulations in Fig. [131 which is to be compared 
with Fig.[TOl The diff'erences are again small and may be mainly 
caused by the temperature diff'erences, higher (average) temper- 
atures leading to fewer cases of multimodal ;^^ distribution. 

As a final deviation from single modified black bodies we ex- 
amined spectra where the spectral index is JS = 2.0 up to 300 yum 
andyS = 1.5 at the longer wavelengths. The break inyS does not 
cause any significant change in the number of double ;^^ minima, 
either in the simulated Planck observations (Fig. [141 vs. Fig- O 
or in the simulated Herschel observations (Fig. [l5lvs. Fig. [TTT) . 



4. Discussion 

The tests with the grey body spectra showed that with four wave- 
lengths 100, 350, 550, and 850 yum and signal-to-noise ratios 
similar to those in Fig. [3 the x^ values of the modified black 
body fits exhibit multiple minima in up to -10% of the cases. 
This is similar to the fraction of 'anomalous' solutions that were 
seen for the radiative transfer models of cylindrical filaments 
(rc<10K)inFigs.[2land[3l 

With the five wavelengths of 100, 160, 250, 350, and 500 yum 
the number of multiple ;^^ minima was lower by a factor of ~ 100 
(see Fig. [TT]) when the 350 yum signal-to-noise ratio was simi- 
lar to that of the previous case with four wavelengths (S/N~20). 
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Fig. 13. As Fig.[T2lbut for simulated Herschel observations. 
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Fig. 14. Frequency of the recognised double x^ minima for the 
modified black bodies with a break in the spectral slope. The 
data correspond to simulated observations at lOOyum and three 
Planck wavelengths. 
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Fig. 15. Frequency of the recognised double x^ minima for the 
modified black bodies with a break in the spectral slope. The 
figure correspond to simulated Herschel data with the 350 yum 
S/N identical to that in the Planck case in Fig.[T4l 



Therefore the better wavelength sampling provided by the set 
of five frequencies is the main reason for the disappearance of 
the double x^ minima. With four wavelengths, the anomalous 
solutions always corresponded to a 350 yum value that was over- 
estimated by ~ 2cr relative to the neighbouring wavelengths (see 
Fig. [6]). With five wavelengths, an anomalous solution might re- 
quire a similar error in two neighbouring channels (e.g., 250yum 
and 350 yum) that, for normally distributed errors and 2 - cr de- 



viations should be less likely by a factor of 50. The number of 
anomalous solutions increased to the 10% level only when the 
S/N ratios were decreased by a factor of six and the 350 yum S/N 
ratio was close to three. 

Are the separate x^ minima of practical importance? In the 
extreme cases like the one shown in Fig. \5\ the multimodality 
of the x^ values causes clear complications. When the minima 
are of similar depth, the result obtained by optimisation meth- 
ods depends on the initial values. Furthermore, an infinitesimal 
change in the flux values or in their weighting can switch the 
global minimum between the low- and the high-temperature so- 
lution. The difl'erence can be several degrees in colour tempera- 
ture and several units in spectral index. Because of the strongly 
non-Gaussian behaviour of the problem the uncertainties de- 
duced from the local shape of the x^ surface will radically un- 
derestimate the true uncertainties of Tq andy^obs- In principle 
the problem can be avoided by calculating the;^^ values over the 
whole parameter plane. This is expensive but may also not yet 
be a completely satisfactory solution. The examples have shown 
that the shape of the x'^ surface can change rapidly as, for ex- 
ample, the error estimates are changed. In Sect. I3.2l we noticed 
that double minima could appear and again disappear when the 
factor ^100 was changed at a 10% level. This means that also the 
uncertainty of the flux uncertainties and their impact on the x^ 
surface should be examined. However, the real error estimates 
are rarely known to a precision of 10%. On the positive side, 
many anomalous solutions may be recognisable from the SED 
plots. In Fig. [5] the low-temperature solution would certainly be 
treated with some caution, in spite of it corresponding to the 
lowest ;^^ value. The situation becomes less clear when the dis- 
tance between the minima is shorter. One should always take 
into account that for low signal-to-noise ratio data there is a non- 
negligible probability, possibly even in excess of -10%, that the 
deepest x^ minimum is not the minimum closest to the true so- 
lution. 

The problem is less tractable in large surveys of individ- 
ual sources or maps with thousands of pixels. It becomes dif- 
ficult to examine each SED fit by eye and the calculation of 
the full x^ surfaces may become impractical. If not accounted 
for, even a few anomalous solutions can significantly aff'ect the 
deduced shape of the I3(T) relation. One can use Monte Carlo 
methods to estimate the number and the influence of the out- 
liers, incl uding the efl'ects of the non-Gaussian statistics (e.g. 
iPlanck Collaboration et al. 201 1 )). Conversely, one can try to 
directly recover the true J3(T) relation with Bayesian methods 
(IVeneziani et al .1120 id: [Kelly et al.ll201 ih. In those cases the mul- 
timodal;^^ distribution may not be a significant complication be- 
cause the methods are aware of the shape of the likelihood func- 
tion and that is additionally modified by the prior. However, an 
accurate knowledge of the statistics of the observed intensities is 
still needed. 

To quantify the role of the double x^ minima (as opposed 
to the general influence of noise) we examined again the simu- 
lated observations of Fig. [51 For each combination of S (350yum) 
and ^100, we took 200 samples of (T^c, y^obs) corresponding to 
the difl'erent temperatures used in the simulation (8.0, 10.0, and 
12.0 K) and a fixed input spectral index of yS = 2.0. The resulting 
linear correlation coefficients between Tq and/^obs are shown in 
Fig.[T6b. In the absence of noise the correlation should be zero. 
The observed correlation varies from -0.5 (for the data with the 
highest signal-to-noise ratio) to — 0.8. Note that the correlation 
coeflftcient does not increase monotonously to the lowest S/N 
ratios because of the strong non-linearity of the ySobs(^c) rela- 
tion. The frame b of Fig. [16] shows the correlation coeflftcients 
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Fig. 16. Linear correlation coefficients between Tq and y^obs 
for the models of Fig. [9] as functions of the source intensity 
S(350fim) and the weight given in the fit to the 100 yum point 
(lower value of ^loo implies a larger weight). Frame b is the 
same after removing the points that correspond to the presence 
of multiple local ;^^ minima. 



excluding those cases where the x^ shows two minima (for that 
particular value of Kiqq). The effect of the cases of double x'^ 
minima is visible at a 10% level. 

The anomalous solutions may be identified by the unre- 
alistic values of the colour temperature and the spectral in- 
dex. However, the presence of multiple x^ minima also directly 
serves as a warning sign. Figure [T7] displays the distribution of 
the (Tq, y^obs) values for the modified black body model with 
r = 12K andy5=1.5 that had particularly many multiple min- 
ima (see Fig. O. The points plotted in Fig. [T7] corresponds to 
the solution obtained from optimisation with initial values close 
to the true solution. The figure shows 500 points for each value 
of S(350/dm), all corresponding to the default weighting with 
^100=1.0. The blue points are the cases where two x^ minima 
were detected with any of the tested ^loo values (see Sect. 13.2b . 
The red points are a subset where double minima were detected 
with the present value of Kioo=l.O. The blue points are seen to 
avoid the locus of the correct solution that still contains most of 
the data. Especially the red points are concentrated in the low- 
temperature tail. Of all the points below the colour temperature 
of 8 K, 42% corresponded to cases with multiple ;^^ minima (or 
extremely shallow single minima that for numerical reasons lead 
to different parameter estimates). If the test is carried out using 
all ^100 factors between 0.5 and 1.5, the percentage increases to 
70%. This suggests that similar tests should be useful also for 
real observations. 



5. Conclusions 

We have studied the behaviour of x'^ fits of SEDs that are ei- 
ther sums of modified black bodies or are based on the ra- 
diative transfer modelling of dust emission from cylindrical 
clouds. Using combinations of wavelengths relevant for the cur- 
rent Planck and Herschel satellite studies, we have examined the 
effect of noise on the shape of the confidence regions in the (Tq, 
ySobs) plane. 

The results have lead to the following conclusions: 

- In addition to the usual symmetric error banana, the x^ dis- 
tribution can exhibit asymmetries of varying strength. The 
expectation values are close to the result obtained in the ab- 
sence of noise, but are not without bias. 

- For low signal-to-noise data (S/N below 10) in the lOOyum, 
350 yum, 550 yum, and 850 yum bands, the noise distribution of 
(rcy^obs) values develops an asymmetric tail that can extend 
to low temperatures and very high spectral indices. 
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Fig. 17. (Tq, JSohs) values for simulated modified black bodies 
with r = 12 K andyS =1.5 (see Sect. 13.2b . The data points corre- 
spond to the normal weighting of data (^loo^l.O) and to all S/N 
ratios shown in Fig.O. The red points denote cases where dou- 
ble ;^^ minima were detected with Kioo=l.O and the blue points 
the other cases where double minima were seen with any value 
of ^100 between 0.5 and 1.5 



- Under the same conditions, the x^ distribution of individual 
measurements can exhibit two distinct minima. A very small 
change in the weighting of the frequency points or in the 
noise can shift the best solution from one minimum to the 
other. This can correspond to a change of several degrees in 
the colour temperature and a change of several units in the 
spectral index. 

- Herschel observations were simulated using five wave- 
lengths, 100, 160, 250, 350, and 500 yum. For the radiative 
transfer models the error distributions remained relatively 
symmetric and very few cases with multiple ;^^ minima were 
detected. This although the signal-to-noise ratios were lower 
than in the previously examined four wavelength case. 

- Investigation of pure modified black bodies (plus noise) 
shows that deviations from a single modified black body, 
such as in the case of line-of- sight mixing of temperatures, 
has no significant effect on the appearance of double ;^^ min- 
ima. 

- The main factor behind the double x^ minima is the noise, 
but the susceptibility depends greatly on the set of wave- 
lengths used. Comparing the four and five wavelength cases, 
equal numbers of double minima where seen when the 
signal-to-noise ratio of the latter were lower by a factor of 
six (S/N~3). 

- The asymmetries or the complete split of the error banana 
have implications for dust studies. It can affect the interpre- 
tation of the observations of individual targets and the relia- 
bility of the I3(T) relations derived from low signal-to-noise 
data. The probability distributions of Tq andy^obs can be non- 
Gaussian, strongly non- symmetric, and possibly even multi- 
modal. These features are sensitive to the assumptions of the 
flux uncertainties and this should be taken into account, even 
in the Bayesian analysis. 
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Appendix A: Second example of bimodal;^^ 
distribution 



Figure IAJI shows another example of a bimodal;^^ distribution. 
The measured values at 100, 350, 550, and 850 yum are 0.14, 
5.01, 3.34, and 1.12MJy sr"^ with the original uncertainties of 
0.06, 0.12, 0.12, and O.OSMJy sr'^. In the modified black body 
fit, when the weight of the lOOyum measurement is varied, the 
solution jumps between a low p and a high p solution. As in 
Fig. [3 the 'correct' solution (i.e., the one closer to the original 
noiseless spectrum) is obtained when the weight of the lOOyum 
point is increased. It is important to notice that in this noise real- 
isation, compared to the noiseless spectrum, the lOOyum is again 
very close to the true value, while the 550 //m is high by more 
than 2(T. 



Appendix B: Tlie bias and scatter of the Tq and ySobs 
values 
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Fig. B.l. Distribution of Tq (red squares and left-hand axis) and 
y^obs (blue circles and right-hand axis) for single modified black 
body spectra with noise and data at the wavelengths of 100, 350, 
550, and 850yum (cf. Fig.O. The solid vertical lines connect the 
quartile points and the open circles corresponds to the 10% and 
90% percentage points of the parameter distributions. 
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Fig. B.2. As Fig. lB.ll but for observations at the wavelengths of 
100, 250, 350, and 500 yum (cf. ¥ig.\M- 



In Sect. 13. 21 we examined the appearance of multiple ;^^ minima 
for modified black body spectra. The parameters Tq andy^obs can, 
of course, have bias and scatter also independent of the problem 
of possible double x^ minima. Figures IB. II and IB. 21 show these 
for the tests with single modified black body spectra. The biases 
and the standard deviations of the colour temperature and the 
spectral index showed no significant dependence on the param- 
eter ^100 and therefore the figures are drawn using only values 
^100=1.0. Otherwise the figures correspond to the cases shown 
inFigs.[9]and[IOl 
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Fig. A.l. Another example of a case with the;^^ having two distinct minima in the (T^c, JS) plane. The jump between the two minima 
takes place between the cases where the original lOOyum error estimate is scaled by 0.702 (left frames) and 0.700 (right-hand 
frames). 
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